Ab initio study of electron transport in dry poly(G)-poly(C) A-DNA strands 
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The bias-dependent transport properties of short poly(G)-poly(C) A-DNA strands attached to 
Au electrodes are investigated with first principles electronic transport methods. By using the non- 
equilibrium Green's function approach combined with self-interaction corrected density functional 
theory, we calculate the fully self-consistent coherent I-V curve of various double- strand polymeric 
DNA fragments. We show that electronic wave-function localization, induced either by the native 
electrical dipole and/or by the electrostatic disorder originating from the first few water solvation 
layers, drastically suppresses the magnitude of the elastic conductance of A-DNA oligonucleotides. 
We then argue that electron transport through DNA is the result of sequence-specific short-range 
tunneling across a few bases combined with general diffusive/inelastic processes. 

PACS numbers: 



I. INTRODUCTION 

Charge transport in DNA has been the subject of in- 
tense research over the past decade [T. This was mainly 
motivated by potential applications in a variety of tech- 
nological fields such as nanoelectronics [21 [3], gene se- 
quencing [4| |5] and chemical sensing [6l. Furthermore 
additional interest has stemmed from the still unclear 
role that DNA electron conduction may play in some bi- 
ological functionalities such as DNA repair [7 . Numerous 
experimental studies on DNA conductivity have yielded 
a range of conflicting results with insulating [8| , semicon- 
ducting [9l [To] , metallic [TTl [12] and even superconduct- 
ing [13] behaviours all being observed. It has been then 
rationalized [M] [15] that the large differences in conduc- 
tivity observed in different experiments may be due to 
the strong dependence of the charge transport process 
in DNA on various intrinsic and extrinsic factors. These 
include DNA conformation, base pair sequence, geomet- 
ric fluctuations, interaction with the solvent and details 
of the experimental setup such as the molecule-electrode 
junction geometry. Nevertheless some consensus is be- 
ginning to emerge that charge transfer in DNA is due 
to a combination of both relatively short-range coherent 
tunneling and longer range incoherent diffusion. 

On the theoretical front different approaches have been 
employed to investigate this complex phenomenon and 
a number of comprehensive reviews may be found in 
the literature [Mj [161 [17] . ^b initio methods are capa- 
ble of providing an accurate description of the detailed 
electronic structure of DNA while treating explicitly sol- 
vent and counter ion interactions [T8H2T] as well as the 
molecule-electrode contact geometry in devices. Classi- 
cal and quantum molecular dynamics simulations have 
yielded much insight into the changes in DNA electronic 
structure induced by structural fluctuations and solvent 
dynamics [20ii21j. However, ab initio transport calcula- 
tions on DNA are currently limited to the linear response 
limit in the static regime because of the large system size 
and the absence of a tractable framework for treating dy- 



namical effects in electron transport from first principles. 

Complementary to ab initio schemes are model Hamil- 
tonian approaches. These are well suited and have been 
employed extensively to investigating dynamical effects 
[22H3QJ such as structural fluctuations and disorder in- 
duced by interaction with solvent molecules, both of 
which are considered crucial to charge transfer in DNA. 
Therefore ab initio and model Hamiltonian investigations 
presently complement each other with the parameters 
for model Hamiltonians being extracted from static first 
principles calculations and efforts to futher develop hy- 
brid approaches underway ^O] Hj . 

Even in the static limit, the full potential of the ab ini- 
tio approach to studying quantum transport in DNA has 
so far not been realized. Most first principles transport 
calculations to date have been restricted to ideal DNA pe- 
riodic structures [32 and studies that explicitly consider 
attachment of DNA oligomers to realistic electrodes have 
been relatively few [33 . Furthermore, transport calcula- 
tions have thus far been carried out only in the zero-bias 
limit. Interestingly, over the past five years a number of 
new experiments, which can be directly compared to ab 
initio calculations, have emerged. These are scanning 
tunneling spectroscopy (STS) measurements on single 
DNA oligomers deposited on Au electrodes and provide 
a systematic mapping of the electronic density of states 
(DOS) around the Fermi energy. Xu et al. |34] employed 
ultrahigh vacuum STS to investigate the electronic prop- 
erties of single thiolated 12-base-pair poly(GC)-poly(GC) 
DNA molecules on a Au(lll) surface at room tempera- 
ture. Shapir et al. [35 used transverse STS measure- 
ments, to map the energy spectrum of poly(G)-poly(C) 
DNA molecules deposited on gold. Both experiments 
produced semi-conductor-like I-V characteristics with a 
well defined voltage gap. 

In this work, we investigate longitudinal (along the 
DNA helix axis) electron transport through single 
poly(G)-poly(C) (pGpC) A-DNA molecules attached to 
Au(lll) electrodes first in the dry and then in moder- 
ately hydrated conditions. The homogeneous 7r-stacking 



and relatively small band gap in pGpC DNA in principle 
provides the best possible conditions for band-like trans- 
port. However, already in dry conditions the electronic 
wave-functions of short pGpC DNA strands (around of 
10 base-pairs in length) are strongly localized by an elec- 
tric dipole along the helical axis. Thus the longitudinal 
transport in dry DNA is dominated by coherent tunnel- 
ing through extremely localized states located at various 
positions along the molecule and the currents are typi- 
cally tiny. For instance we calculate the entire I-V curve 
of a 6 base-pair pGpC DNA oligomer attached in var- 
ious ways to the electrodes and show that the currents 
originating from coherent electron tunneling at voltages 
of around 1 Volt are several orders of magnitude smaller 
than those observed experimentally through single DNA 
transport devices. 

Wetting the DNA strand does not improve the coher- 
ent transport picture. Water and counter- ions in fact 
produce a dipolar field, which partially screens the na- 
tive DNA dipole. This generates additional electrostatic 
disorder. The DNA molecular orbit als responsible for 
the electron transport remain largely localized but now 
their actual spatial distribution and energy position with 
respect to the electrodes' Fermi level depends on the par- 
ticular geometrical configuration of the solvent. The as- 
sociated transmission coefficient as a function on energy 
is then qualitatively similar to that of the dry case. How- 
ever, the presence of the solution produces an inhomoge- 
neous level broadening that we calculate to be as large 
as 1.5 eV and provides an efficient channel for inelastic 
transport. From this analysis emerges a picture of elec- 
tron conduction in DNA dominated by coherent electron 
tunneling on a length-scale of a few base-pairs combined 
with incoherent diffusion regulated by the solvent. 

The paper is organized as follows. In the next two sec- 
tions we will introduce our theoretical methods and dis- 
cuss the basic electronic structure characteristic of pGpC 
DNA oligomers both in isolation and when attached to 
Au electrodes. This will set the most appropriate elec- 
tronic structure theory for the problem. In particular we 
will show that correcting the local density approximation 
(LDA) for electron self-interaction [37 is fundamental in 
this case in order to obtain a realistic molecular level 
alignment. Then we will move to the transport calcula- 
tions by looking first at the zero-bias limit and then to 
finite bias situations. Finally we will consider the effects 
of solvation and conclude. 



Hartree-Fock total energy calculations. When Siesta 
is employed, the LDA with the Ceperly-Alder parame- 
terization for the correlation energy [40] is used as the 
exchange-correlation functional. Self interaction correc- 
tions over the LDA are incorporated within the atomic 
self interaction scheme (ASIC) [41]. Standard norm- 
conserving scalar relativistic pseudopotentials with the 
following reference configurations are used for the differ- 
ent atomic species: Rls^.C 2s'^2p^ , N 25^2^^, O 25^2^'^, 
Na 35\ P 35^3^^, S Zs'^Zp^, and Au 65^ An optimized 
double-C basis set is employed for all the atomic species 
except for Au where a single-^ basis is used. Sampling in 
real space is done through a uniform grid with an equiv- 
alent cutoff energy of 150 Ry. When using Gamess the 
SBKJC effective core potentials of Stevens et al. [42 are 
employed together with the associated SBKJC basis set. 
GGA and hybrid-DFT calculations are carried out using 
the PEE [43^ and PBEO [44] exchange correlation func- 
tionals respectively. 




FIG. 1: (Color online). Geometric configuration of A- 
DNA both in its free form and when attached to the elec- 
trodes: (a,b) front- and side-view respectively of a 11 base- 
pair poly(G)-poly(C) A-DNA double helix, (c) Two-probe 
device set up used for transport calculations where in the 
DNA molecule is attached to Au electrodes. Color code: 
Yellow=Au, Green=C, Red=0, Orange=P, Blue-Gray=N, 
Cyan=H, Pale-Green=S. 



II. METHODS 

The Siesta [38^ code, which employs a numerical lo- 
calized orbital basis set and norm-conserving pseudopo- 
tentials forms our primary computational platform. In 
fact most of the density functional total energy calcula- 
tions and geometry optimizations are performed within 
Siesta. In addition the Gamess [39] quantum chemistry 
package is employed to carry out some hybrid-DFT and 



Transport calculations are performed with the 
Smeagol code [45J [46], which combines the non- 
equilibrium Green's function method (NEGF) J47l l48]. 
with DFT as implemented in Siesta [38]. The various 
simulation parameters listed above for Siesta are then 
transferred to the transport calculations. The scattering 
region of the transport simulation setup comprises, in 
addition to the molecules, five atomic layers of Au(lll) 
included at both sides of the DNA. These are sufficient to 



adequately screen charging at the molecule- Au interface. 
In order to reduce the system size the Au-5d shell is left 
in the core and only the 6s electrons are treated as va- 
lence. We use an optimized single-^ basis for the Au-6s 
orbitals, specifically tuned to reproduce the correct work 
function of the Au(lll) surface. This is an established 
and well documented procedure already used with suc- 
cess in the past for describing electron transport across 
small molecules [49] , 

The typical size of our simulations and the intrin- 
sic weak coupling between the DNA molecules and the 
electrodes demand particular care when integrating the 
Green function yielding the charge density. Such an inte- 
gral is split into two parts, one over the complex energy 
plane and one along the real axis [4F . The complex inte- 
gral is performed over a uniform mesh of 512 imaginary 
energies, while for the real part we have implemented a 
new mesh refinement algorithm, necessary to integrate 
the extremely sharp features of the DOS [50 . Typically, 
in order to integrate the non-equilibrium density over a 
2 eV window, the final mesh includes 8000-10000 energy 
points with the denser energy spacing being ^ 10~^ eV 
around the sharpest transmission peaks. Furthermore, 
in order to boost the convergence for large system sizes, 
a number of technical solutions have been implemented 
including a singularity removal procedure to evaluate ac- 
curately the leads self-energies [51^ . 



III. DFT ELECTRONIC STRUCTURE 
A. Isolated molecules 

The system studied first is an eleven base-pair pGpC 
DNA molecule with the A-DNA conformation. This con- 
stitutes one complete twist of the DNA double helix 
structure (see figure fl]), i.e. in solid state terms it forms 
the A-DNA one-dimensional unit cell. The molecule, 
consisting of 715 atoms, is 30.5 A long and has an aver- 
age twist angle of 32.72^. The negatively charged DNA 
Sugar-Phosphate backbone is neutralized by adding two 
H ions per base-pair to the PO4 groups along the length 
of the backbone. By imposing periodic boundary condi- 
tions on this system, an infinitely long DNA chain can be 
simulated. Alternatively, by suitably pacifying the chain 
at both ends and by relaxing the positions of the atoms 
at the edges, a finite 11 base-pair oligomer may be mod- 
elled. For the periodic structure, we sample the Brilluion 
zone at two /c-points but even a F-point calculation would 
be adequate. We present the electronic structure of both 
the periodic and finite chains in order to investigate the 
differences between the two. 

Figure |2] shows the calculated DOS projectd onto dif- 
ferent parts of the DNA chain and displays separately 
the contributions from Guanine, Cytosine and the Sugar- 
Phosphate (SP) backbone groups corresponding respec- 
tively to the Guanine (G-SP) and Cytosine (C-SP) side 
of the helix. In the periodic chain (left panel of figure ^ 
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FIG. 2: (Color online) Density of states projected onto differ- 
ent constituent parts of poly(G)-poly(C) A-DNA chains: (left 
panel) projected DOS for the infinite periodic DNA structure 
with 11 base-pairs in the unit cell; (right panel) projected 
DOS for a finite 11 base-pair oligomer. The purple vertical 
line indicates the Fermi level. 



a well defined band-like DOS appears and the individual 
bands are well separated from each other. The highest 
occupied (HO) band is derived from Guanine and has 
a rather narrow bandwidth of about 60 meV, while the 
lowest unoccupied (LU) band corresponds to a set of Cy- 
tosine states with a relatively larger bandwidth of 300 
meV. The system has an overall (LDA) bandgap of 2 eV. 
Among the occupied bands, the first three highest ones 
are derived from orbitals localized on Gunaine and there 
is little Cytosine or backbone related DOS upto 1.5 eV 
below the HO band. Among the empty states, the first 
Cytosine-derived band is followed by another originat- 
ing from orbitals localized at the Sugar-Phoshate groups 
attached to the Cytosine chain. 

When going from the periodic to the finite molecule 
(right panel of figure [2|, each band splits up into 11 
discrete levels and the derived energy eigenvalues now 
spread across a much wider energy range. For instance, 
individual Guanine-derived eigenstates are spread over a 
400 meV range on both sides of the 60 meV HO band 
previously calculated for the periodic case. As a result, 
the highest occupied molecular orbital (HOMO) and the 
lowest unoccupied molecular orbital (LUMO) are now 
separated by ~1.25 eV. In order to explore the effects of 
having finite fragments instead of infinite periodic DNA 
chains on the electronic density we plot the charge den- 
sity corresponding to different representative eigenstates 
(see figure [3|. On the one hand the eigenstates from the 
HO and LU bands of the periodic system are seen to 
spread uniformly across the base pair stacking along the 
Gunaine and Cytosine chains respectively. On the other 
hand, the eigenstates of the finite molecule are seen to 
be spatially localized on one or two bases at most. 

This localization of the wavefunctions is the result of 
the charge dipole along the helical axis created when 
the periodic DNA chain is truncated to form the finite 
molecule. The dipole points from the 3' end towards the 
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FIG. 3: (Color online). Projected DOS (top panels) and charge density distribution (lower panels) for both infinite (left) and 
finite (right) DNA strands. The projected DOS are displayed over an energy window showing the first Guanine and Cytosine 
levels around the Fermi level. The arrows and the corresponding numbers mark the various energy levels, whose corresponding 
charge density distributions are displayed in the lower panels. Note the delocalized nature of the HO and LU bands of the 
infinite chain in comparison with the localized states of the finite fragments. 



b' one of the Guanine chain and it Stark-splits the energy 
levels corresponding to the various bases. In fact, for the 
infinite chain the average electrostatic potential at each 
individual base is the same for all the bases leading to the 
formation of extended states albeit with a small band- 
width since the interbase hopping is weak. In contrast 
for the finite chain, the dipole shifts up (down) the av- 
erage electrostatic potential of the Guanines (Cytosines) 
at the 5' end with respect to those at the the 3^ Since 
such a potential shift is larger than the effective hopping 
between the bases, the electronic wavefunctions become 
localized. The HOMO is now spatially placed at the 5' 
end of the Guanine chain and states with progressively 
lower energies are arranged in order from the 5' end to- 
wards the 3' one. Similarly, the LUMO is now localized 
at the 5' end of the Cytosine chain. This type of dipole- 
induced localization is expected to be more pronounced 
at the edges and less in the middle of the molecule es- 
pecially as the chain length increases. Thus short pGpC 
DNA fragments in dry conditions, even in the absence 
of other forms of static disorder, support only localized 
wavefunctions because of the electrostatic dipole along 
the helical axis. 



B. DNA fragments attached to Au electrodes 

We now describe the scattering region setup for a 11 
base-pair pGpC strand. DNA strands are connected to 
the Au(lll) surface via standard (CH2)3-SH thiol link- 
ers. In the reference geometry, the thiol linkers are at- 
tached at both the 3' and 5' ends of the Guanine and 
Cytosine chains providing a total of 4 linkages between 
the molecule and the Au surface. We label this geometry 
G3'C3'G5'C5'. Connection to the Au(lll) electrodes is 
made as follows. First, the S atom of one of the linker 
groups is placed at a hollow site forming an optimal bond 



length to the surface Au atoms [53 . The molecule is then 
re-oriented so that the distance from the Au surface of 
the S atom on the second linker group is optimized. This 
procedure fixes the orientation of the molecule. The sec- 
ond Au electrode is then rigidly shifted to optimize its 
distance from the two thiol linker groups on the opposite 
side of the molecule. The device simulation cell includes 
five Au atomic planes on each side of the molecule. The 
lateral dimensions of these are those of a 9 x 7 supercell 
constructed from a rectangular unit cell along the fee 
(111) direction. The unit cell for the electrodes is also 
a 9 X 7 supercell. A lattice constant of a = 4.1812 A is 
used for Au. Periodic boundary conditions are enforced 
in the plane transverse to the transport direction. The 
final simulation cell contains 2017 atoms, of which 1260 
are Au atoms of the electrodes. The basis set expands 
over 5320 local orbitals. 

Next we look at the electronic structure of DNA 
strands attached to Au electrods. In NEGF transport 
calculations it is fundamental to describe accurately the 
ionization potential (IP) and electron affinity (EA) of 
the molecule as well as the work function, W, of the elec- 
trodes. This allows one to reproduce the correct align- 
ment of the molecular levels with respect to the elec- 
trodes' Fermi level, E^. In the case of pGpC this task 
is complicated by the fact there is no experimental infor- 
mation about the IP and EA of pGpC fragments. Our 
strategy is then that of employing a DFT functional that 
reproduces well the experimental IP of Guanine and Cy- 
tosine in their gas phase and then to calculate that of 
pCpG fragments as a function of the number of base- 
pairs in the fragment. The parameter- free PBEO hybrid- 
DFT functional has been shown previously to reproduce 
molecular excitation energies in good agreement with ex- 
periments [44]. For Guanine (Cytosine) the calculated 
IPcai = 8.11 eV (IPcai = 8-65 eV) compares well with 
the experimental one IPexp = 8.2 eV (IPexp = 8-9 eV) 




12 3 4 5 6 7 
Number of base-pairs (n) 



-10 1 

Energy (eV) 



1 

Energy (eV) 



FIG. 4: (Color online). DNA electronic structure calculated at various levels of approximation and level alignment between 
the molecule and the Au electrodes. In the left panel we show the negative of the ionization potential (solid lines) and of the 
electron affinities (dashed lines) calculated from PBEO-ASCF and Hartree-Fock for pGpC DNA strands of different lengths. 
Also included for comparison are the highest occupied and lowest unoccupied eigenvalues from PBE calculations and the Au 
workfunction. In the middle panel we present the atom projected DOS obtained respectively from LDA and ASIC calculations 
for the 11 base-pair DNA molecules attached to Au electrodes with the reference G3^C3^G5^C5^ geometry. Finally the right 
panel displays the zero-bias transmission coefficients as calculated with LDA and ASIC for the same system. The Au Fermi 
level is set to eV. 



[5F. Here both the IP and EA are calculated with the 
ASCF scheme. The vertical IP of the molecule in vac- 
uum is computed as IP=E{N — 1) — E{N)^ which is the 
energy difference between the singly positively charged 
and the neutral molecule {N is the total number of elec- 
trons at neutrality). Similarly, the EA is calculated as 
EA=E{N) - E{N + 1) which is the energy difference 
between the neutral and the singly negatively charged 
molecule. 

In the first panel of figure [4J we plot the ASCF- 
calculated negative IP (green solid line) and the nega- 
tive EA (green dashed line) for DNA strands ranging 
from 1 to 8 base-pairs in length. Interestingly we observe 
that both the quantities present a substantial length- 
dependence. The IP decreases as the number of pairs 
increases and finally converges towards a value of around 
5.65 eV for the 8 base-pair chain. In contrast the EA has 
an opposite behavior, it reduces as the number of base- 
pairs increases and it approaches the value of 1.2 eV for 
8 base-pairs. 

The negatives of the IP and EA set the alignment of 
the molecular eigenvalues with respect to the vacuum 
level. Furthermore, the true HOMO-LUMO gap of the 
molecule is simply EA-IP and it is calculated to be to be 
4.47 eV. In figure [4] we also plot (blue lines) the removal 
and addition energies estimated from Hartree-Fock (HE) 
calculations by using Koopman's theorem (i.e. from the 
position of the single particles eigenvalues) . HE calcula- 
tions usually tend to overestimate the IP and the EA of 
DNA strands. This appears to be the case here with the 
IP converging to around 6.61 eV for the longer chains. 



which is roughly 1 eV higher than the ASCF PBEO esti- 
mate. Similarly, also the EA is underestimated with the 
LUMO eigenvalues coming out positive (i.e. the addi- 
tional electron is not bounded). 

Although the IP and EA can be estimated with good 
accuracy from ASCF calculations the corresponding sin- 
gle particle energy levels, enoMO and eLUMO respectively, 
which are the ones that enter in a NEGE calculation, are 
not readily available in DFT. This is because chomo ob- 
tained from semi-local DFT functionals is typically too 
high in energy i.e. it underestimates the IP. Likewise 
clumo is usually too low. Again this is the situation 
encountered here for the PBE-GGA as demonstrated by 
the results of figure [4] (red lines). Clearly, enoMO from 
PBE is too high in energy compared to -IP estimated 
from ASCF. The work function of the Au(lll) surface 
is calculated to be '^5.3 eV and it is indicated by the 
orange dotted line in figure [4] We then need to cor- 
rect the DFT single particle levels in such a way for the 
chomo to correspond to the negative of the IP as calcu- 
lated with the ASCF scheme. This is achieved by using 
the ASIC method ^, a strategy that has already been 
successfully employed in earlier molecular transport cal- 
culations [491151]. 

In general, the ASIC scheme uses a scaling parame- 
ter, a, ranging between and 1 (a = reduces to the 
LDA, a = 1 corresponds to the full SIC for 1 electron 
systems). The value of a used is related to the charge 
screening properties of the material under consideration 
and it is typically ~1.0 for small molecules, around 0.5 
for insulating oxides and vanishes for metals ^41,. The 



electronic spectra of numerous wide gap semiconductors 
have been studied previously using the ASIC scheme with 
a value of a=0.5 [41 . For the case of the DNA oligomers, 
we choose a value of a that produces an enoMO approxi- 
mately equal to the negative IP from ASCF calculations. 
For this choice an a in the range of 0.5-0.6 is found to be 
appropriate. 

In figure [4] (middle panel), we present the DOS for 
the scattering region in the reference geometry calcu- 
lated using both ASIC {a = 0.6) and LDA (note that 
LDA and PBE eigenvalues are practically the same). In 
LDA enoMO is pinned at the Au £^f- This happens be- 
cause within LDA/GGA the gas phase chomo of the 
DNA strand lies above Ep resulting in the creation of 
a charge dipole at the metal- molecule junction, which 
acts to align chomo with the Au work-function. In fact 
the sharp peak at Ep seen in the LDA DOS is comprised 
of 11 nearly degenerate levels corresponding to the 11 
base pairs. These states, whose energies spread across 
a much wider range when the molecule is isolated, are 
all brought into near resonance because of the fractional 
charging of the molecule. In contrast ASIC places E^ 
in the DNA HOMO-LUMO gap at approximately 0.7 eV 
from enoMO- Since chomo is already below Ep and the 
first empty molecular state eLUMO is well above it, no 
additional charge dipole is created. Also the structure of 
the HOMO manifold is not altered with respect to the 
isolated molecule. 

Note that although enoMO is aligned with the nega- 
tive IP within ASIC, clumo does not correspond to -EA 
which lies higher. This is formally not a problem as eigen- 
values of empty molecular levels need not equal electron 
addition energies within DFT. However, the fact that the 
HOMO-LUMO gap is underestimated with respect to the 
quasiparticle gap within the current approach should be 
borne in mind while interpreting results from transport 
calculations. The picture from the DOS is also carried 
over into the transmission coefficients calculated from 
LDA and ASIC (right panel of figure H). Becaue of the 
spurious pinning of enoMO at E^ LDA predicts a high 
transmission at zero bias while ASIC predicts a semi- 
conducting system with very little transmission around 
E — Ep. 



IV. ZERO BIAS RESULTS 

We next move to analyzing in detail the transport 
properties of pGpC, by starting from the zero-bias limit. 
From here on, all the transport results presented as well 
as the discussion that follows will be based on ASIC cal- 
culations. In figure [5] we present the zero-bias transmis- 
sion coefficient, T{E)^ as a function of energy, E^ for 
the reference geometry (11 base-pairs and connections to 
gold through both the strands). The T{E) is typical of a 
tunneling system and presents rather small transmission 
values in the gap region away from both the HOMO and 
LUMO (note the T{E) is plotted in logarithmic scale). 



Furthermore, the broadening and maximum height of the 
different transmission peaks, corresponding to the differ- 
ent molecular orbit als, vary significantly across several 
orders of magnitude. In general, the coupling strength 
to the electrodes of DNA molecular orbitals that are de- 
rived mainly from the Guanine and Cytosine bases is very 
weak. This is because the bases are well separated from 
the Au surface by the linking thiol groups as well as the 
sugar phosphate backbone. As a result the coupling in- 
duced broadening of the transmission peaks is very small, 
typically less than ~ 10~^ eV for the HO and LU main- 
fold of states. The resonances corresponding to chomo 
and clumo are indicated by arrows and the correspon- 
dence between the DOS and T{E) is clearly evident. 
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FIG. 5: (Color online). Zero-bias transport properties for 
an 11 base-pair pGpC in the reference G3^C3^G5^C5^ geom- 
etry. In the top left panel we present the projected density 
of states. The HOMO and LUMO are indicated by arrows 
and labelled by 1 and 2 respectively. In the bottom left panel 
the zero-bias transmission coefficient as a function of energy 
is presented. In the two right panels we show the charge den- 
sity distribution corresponding to the HOMO (top) and to the 
LUMO (bottom). Note that these are localized respectively 
near the 5^ end of the Guanine stack and near the 5^ end of 
the Cytosine stack. 

Depending upon the strength and symmetry of the 
coupling of the various molecular levels to the electrodes, 
two types of transmission peaks can be identified. Molec- 
ular orbitals that are localized on one side of the scatter- 
ing region couple relatively strongly but asymmetrically 
to the leads. For example, the electronic density for the 
HOMO level (right panel of figure |5| is well localized 
around the 5' end of the Guanine chain near the left 
hand-side electrode. Therefore, the HOMO couples far 
more strongly to the left hand lead than it does to the 
right hand-side one. This asymmetric coupling leads to 
an attenuated transmission peak at chomo [50 . Simi- 
larly the LUMO charge density is localized on the right 
hand-side at the 5' end of the Cytosine chain thus couples 
more strongly to the right hand-side lead. In contrast, 
orbitals that are localized near the middle of the DNA 
chain couple weakly but more symmetrically to the elec- 
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FIG. 6: (Color online). Zero-bias transport properties for three different bonding geometries of a pGpC 11 base-pairs DNA 
strand attached to gold. A schematic depicting the geometry of three diffrent molecule-electrode linking configurations is show 
in the top panels. All the configurations are obtained from the reference (G3^C3^G5^C5^) one, by removing two thiol linkers. 
The corresponding T(E) for each geometry are presented in the lower panels. In all the figures the light brown line is the 
transmission for the reference G3^C3^G5^C5^ geometry. 



trodes and thus display corresponding peak heights that 
are closer to the maximum transmission of 1. As dis- 
cussed earlier, the eleven Guanine derived states in the 
highest occupied manifold that are sequentially lower in 
energy are arranged from left to right along the Guanine 
chain. Orbitals corresponding to energy eigenvalues near 
the center of the HO manifold are localized near the mid- 
dle of the chain. Furthermore, these states are expected 
to be relatively less localized in space as their correspond- 
ing energy eigenvalues are closely packed. This further 
improves the symmetry of their coupling to the leads. A 
similar analysis also applies to the LUMO manifold. 

Having looked at T{E) for the reference G3'C3'G5'C5' 
geometry, we now investigate the effect of changing the 
molecule-electrode linkage geometry. Accordingly, we 
build three different junctions in which thiol linkers from 
the G3'C3'G5'C5' configuration are selectively removed 
and the ends of the DNA chain suitably pacified with H 
atoms. Furthermore, in removing the alkyl-thiol linker 
groups, ghost basis orbitals constructed to match the 
orbitals on the atoms comprising the linker groups are 
substituted in place of the removed atoms. Note that 
in making these new configurations, the position of the 
DNA chain in space is kept identical to that in the refer- 
ence setup. In figure|6) we show schematics of the geome- 
tries as well as the corresponding zero bias transmission 
coefficients. 

In the G3'G5' configuration, the Guanine chain is at- 
tached to Au by thiol groups at both ends, while the 
Cytosine chain is simply pacified and no linkers are con- 
nected. Similarly in the C3'C5' configuration only the 
Cytosine strand is attached to Au by thiol linkers at 
both ends while the Guanine chain is left unattached. 
In both these situations the electrodes are electronically 
connected only through one of the two strands of the dou- 
ble helix. In contrast the G5'C5' configuration is made 



by connecting the 5' ends of both the Guanine and Cyto- 
sine chains to Au via thiol linkages, while leaving the 3' 
ends free. In this case there is no individual strand which 
is connected at both ends to the electrodes, i.e. electrons 
must be transferred across the two strands in order to 
sustain a current. 

The T{E) plots for these three new geometries are pre- 
sented in figure [6J where we also reproduce again the 
zero-bias transmission of the G3'C3'G5'C5' reference ge- 
ometry (light brown line) for ease of comparison. We 
note that in the G3'G5' geometry T{E) for the Cyto- 
sine derived LU manifold of states is suppressed with 
respect the G3'C3'G5'C5' configuration while the trans- 
mission for the Guanine derived HO manifold is virtually 
unchanged. This is different from the C3'C5' junction, 
where the transmission corresponding to the HO set of 
states is seen to be similarly suppressed, while that for 
the LU manifold remains largely unchanged with respect 
to the reference geometry. Finally, for the G5'C5' con- 
figuration the transmission is reduced for both the HO 
and LU manifolds although the largest reduction is for 
the HO manifold. 

These results show firstly that, even though the thiol 
linkers are not bonded directly to the Guanine and Cyto- 
sine bases but are only attached to the sugar-phosphate 
backbone, they play a critical role in the transmission. 
In fact the electron tunneling from the electrodes to 
the bases is much larger across these organic linkers 
than through vacuum. Secondly, electron tunneling from 
molecular orbitals to the electrodes, even for those or- 
bitals that are localized near one end of the DNA strand, 
effectively occurs through the base-pair stacking. Thus 
in the G5'C5' geometry, the broadening of a T{E) peak 
corresponding to an orbital localized at the 5' end is re- 
duced even though only the thiol linker at the 3' end is 
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FIG. 7: (Color online) Projected density of states and charge density plots for two different device setups comprising 6 base-pair 
pGpC DNA strands attached to Au electrodes. In the top (bottom) left panel we show the charge density corresponding to the 
HO orbital on the DNA molecule connected in the C3^G3^ (G5^G3^) geometry. In the central top (bottom) panel the projected 
DOS for the C3^G3^ (G5^G3^) setup is displayed. The HO and LU molecular states are indicated by arrows. Finally in the top 
(bottom) right panel we present the charge density corresponding to the LU orbital on the DNA molecule connected in the 
C3^G3^ (G5^G30 geometry. 



removed. 



V. FINITE BIAS RESULTS 

We now turn our attention to the transport properties 
of pGpC at finite bias. Since finite bias calculations are 
considerably more computationally intensive than those 
at zero-bias, we now consider a shorter DNA molecule, 
namely a 6 base-pair pGpC DNA strand. The construc- 
tion of the scattering region is similar to the one used 
for the zero bias calculations except for the length of 
the molecule. The lateral dimension of the electrodes is 
the same as before and periodic boundary conditions are 
used in the plane. The connection of the molecule to 
the Au electrodes is once again through (CH2)3-SH thiol 
linkers. We construct two different scattering regions dif- 
fering mainly in the way the molecule is attached to the 
electrodes. In the first configuration labelled C3'G3^ the 
3' ends of both the Cytosine and Guanine chains are con- 
nected to Au via thiol linkers while the 5' ends are left 
free. In the second configuration, labelled G5'G3^ both 
the 5' and 3' ends of the Guanine strand are connected to 
Au while neither end of the Cytosine chain is connected. 
Both these configurations are shown in figure [7J 

The atom projected DOS (atom PDOS) for the two 
systems calculated within ASIC {a=0.6) is also shown in 
Fig.[7| where both the HOMO and LUMO are marked by 
arrows. Additionally, the electron density corresponding 
to the HOMO and LUMO is also plotted. The PDOS 
is qualitatively similar to that of the 11 base-pair case 
with the nature and relative energy ordering of the dif- 
ferent molecular orbitals remaining the same. As before. 



the HOMO charge density is localized at the 5' end of 
the Guanine chain while the LUMO charge density is 
localized at the 5' end of the Cytosine chain. Thus the 
strength and the level of symmetry of the coupling to the 
electrods of different orbitals are expected to vary along 
the base-pair stacks. 

The calculated I-V curve for the C3'G3' geometry is 
shown in figure |8] Note that the absolute value of the 
current is plotted in a logarithmic scale. The magnitude 
of the current is seen to be very small reaching up values 
in the range of 10~^^-10~^^ Amperes for bias voltages be- 
tween 1.2 and 1.4 Volt. These currents are several orders 
of magnitude smaller than those measured in transport 
experiments on single DNA oligomers at comparable volt- 
ages [9l [Ml ES] . The potential drop across the device at 
two different applied voltages is also shown from which 
it is apparent that the voltage drop occurs almost en- 
tirely across the molecule and is flat within the Au layers 
included in the scattering region. The bias dependent 
transmission coefficient is also plotted in figure |8] for three 
different bias points. At zero bias the leads' Fermi level 
lies in the HOMO-LUMO gap of the molecule, where 
T{E) is extremely small. For positive bias voltages the 
Fermi level of the left lead (E^) is raised while that of the 
right lead (^p) is lowered in energy. For negative bias 
voltages the two Fermi levels are moved in the opposite 
direction. Since the native dipole on the 6 base-pair DNA 
strand points from the 3' end towards the 5' end of the 
Guanine chain, when the 5' end is connected to the left 
electrode, the external electric field for positive bias acts 
in the same direction as the native dipole. For negative 
bias on the other hand the external field opposes the na- 
tive dipole. As a result, the the transmission peaks for 
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FIG. 8: (Color online) Finite bias transport properties of the C3^G3^ molecule-electrode configuration. In the left panel we 
show the I-V curve. Note that the current is displayed in logarithmic scale and that for 1/ < we actually show the negative 
of the real current (— /). In the central and right panels we present respectively the transmission coefficients as a function of 
energy at different bias voltages and the potential drop across the scattering region. In the T(E) plots the vertical dashed lines 
mark the bias windows. 
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FIG. 9: (Color online) Finite bias transport properties of the G5^G3^ molecule-electrode configuration. In the left panel we 
show the I-V curve. Note that the current is displayed in logarithmic scale and that for y < we actually show the negative 
of the real current (— /). In the central and right panel we present respectively the transmission coefficients as a function of 
energy at different bias voltages and the potential drop across the scattering region. In the T(E) plots the vertical dashed lines 
mark the bias windows. 



the HO manifold spread out over a larger energy range 
for positive bias while the same set of states are pushed 
closer together into near resonance for negative bias. The 
LU manifold of states also show similar behaviour. Fur- 
thermore, for both positive and negative bias, the trans- 
mission peak corresponding to the HO level is the first 
near to the bias window and therefore the voltage gap in 
this case does not depend upon the position of the LU 
level. At 1.2 Volt for positive bias for example the trans- 
mission peak corresponding to the HOMO is inside the 
bias window while that of the LUMO is still about 0.4 
eV away from it. Similarly for negative bias, the HOMO 
is within 0.1 eV of the bias window while the LUMO is 
much farther away. The overall size of the voltage gap 
for the C3'G3' system is seen to be approximately 2 Volt. 
Such a gap is in good agreement with experimental data. 

In figure |9J we show the I-V curve, T{E) for different 
bias voltages and the potential drop across the scatter- 
ing region for the G5'G3' configuration. Qualitatively, 
the results for the G5'G3' system are similar to those of 
the C3'G3' one. The current is again very small, in the 
range of 10~^^-10~^^ Amperes for reasonable bias volt- 
ages around 1.0-1.5 Volt. Moreover, it is the HOMO level 
that once again conducts first and sets the size of the volt- 



age gap at about 2.2 eV. Thus in spite of marked differ- 
ences in the molecule-electrode connection geometry, the 
transport properties of the two device configurations are 
similar. This is because for the range of voltages inves- 
tigated here the nature of the I-V is determined mainly 
by a feature common to both geometries, i.e. by the fact 
that the HOMO is localized at one end of the Guanine 
chain and it is more strongly coupled to the left lead. 

The extremely small magnitude of the calculated cur- 
rents at sizeable voltages suggests that even in relatively 
short DNA strands, the coherent tunneling mechanism is 
very weak and dynamic or inelastic processes must there- 
fore dominate the charge transfer process. Additionally, 
the large contact resistance across the (CH2)3-SH linkers 
in our simulation geometry might also contribute to sup- 
press the current in comparison to experiment. Although 
in-principle (CH2)3-SH linkers are standard in experi- 
mental setups m [34l [35] , one cannot rule out potential 
imperfections in the junction geometries that allow some 
parts of the DNA strands to make direct contact with the 
electrodes thus reducing the contact resistance. Because 
of their substantial length and large HOMO-LUMO gap, 
the contact resistance across the (CH2)3-SH linkers is 
expected to be high and indeed the effective coupling to 
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the electrodes of even the DNA bases at the ends of the 
strand is seen to be rather weak. From the fuh-width at 
half- maximum, 7HOMO7 of the transmission peak corre- 
sponding to the HOMO orbital localized at the 5' end 
of the Guanine stack, it is possible to make an estima- 
tion of the saturation elastic current [50 for a junction 
with just one base-pair directly attached to Au electrodes 
across (CH2)3-SH linkers, assuming that the coupling to 
both electrodes is equal to 7homo- We find 7homo to be 
^3xlO~^eV which yields a value of ~0.3 nA for the sat- 
uration elastic current [50 through such a junction with 
just one base pair contacted via (CH2)3-SH linkers. Thus 
solvent mediated in-elastic effects might already play a 
role in the process of electron hopping between the elec- 
trodes and the terminal bases of DNA strands. 



VI. SOLVENT EFFECTS 

It is well known that the conductivity of DNA depends 
strongly on environmental effects, in particular on the 
presence of a solvent and additionally counter-ions. The 
solvation state can infiuence the electronic structure of 
DNA in several ways. For instance the confirmation as- 
sumed by DNA double helices is itself determined by the 
extent to which the molecules are surrounded by water 
and by the energy associated with hydration. Addition- 
ally, as shown recently by Rungger et at. [56 , in molecu- 
lar junctions, solvation by water effectively produces elec- 
trostatic gating which can change the average alignment 
of the molecular energy levels with respect to the leads 
E^. Finally, around room temperature the thermal mo- 
tion of water molecules can produce rapidly fiuctuating 
local dipolar fields that in turn lead to fiuctuations in the 
energy level positions of individual molecular orbitals on 
the DNA. Such fiuctuation may play a crucial role in 
charge transport processes through DNA chains [57] , 

In most STS-based experimental studies on charge 
transport through DNA [34, 35 the double stranded 
oligonucleotides are assembled in solution before under- 
going a drying process aimed at removing the solvent 
molecules. However, the exact solvation state of the DNA 
molecules resulting from such a process is usually un- 
known and it is likely that a thin solvation layer still 
surrounds the DNA with the polar water molecules in 
the solvent hydrating both the negatively charged phos- 
phate backbone of the DNA as well as any counter-ions 
that may be present. Therefore, in the context of study- 
ing electron transport through A-DNA attached to gold 
electrodes, we also investigate the effects of a thin layer of 
solvent consisting of water molecules and Na+ counter- 
ions. We note that nominally wet DNA in the B and 
Z conformations has been studied previously by several 
theoretical groups [HI [20l [33l [36] both within a static pic- 
ture as well as including the dynamics of water molecules. 
In particular, Mantz et. al [36 studied the interplay 
between solvent geometry and hole charge localization 
in Adenine-Thymine bridges using QM-MM simulations 




FIG. 10: (Color online) A 6 base-pair pGpC molecule in wet 
conditions. In (a) we show the scattering region setup. This 
consists of a 6 base-pair pGpC A-DNA strand surrounded by 
a few layers of solvent including water molecules and Na^ 
counterions. In (b) we show the planar average of the elec- 
trostatic potential across the system at zero-bias. The green 
line shows the electrostatic potential in the absence of the 
solvent while the black line shows the potential averaged over 
100 time snapshots of the system taken from the classical MD 
simulations. The gray lines show the electrostatic potential 
for the individual snapshots. 



demonstrating the gating effect due to solvent rearrange- 
ments on the pico-second time scale. 

Here we are mainly interested in the differences in the 
electronic structure between the dry and solvated A-DNA 
strands as calculated within the ASIC method. In partic- 
ular we focus on the fiuctuations in the electronic energy 
level positions induced by interactions with the polar sol- 
vent. The simulation setup is shown in figure [Tol We 
consider a pGpC 6 base pair A-DNA duplex attached to 
Au electrodes in the C3'G3' geometry as described earlier 
but now surrounded by a few layers of water molecules in 
the presence of Na+ counter-ions. The strategy adopted 
for investigating the effects of the solution on the elec- 
tronic transport is the same reported recently in reference 
[56] . It consists first in performing classical molecular dy- 
namics (MD) simulations in order to describe the average 
geometrical arrangement of the solvent around the DNA, 
and then in evaluating the average transmission over a 
number of non-time-correlated configuration snapshots. 

In the dry condition the DNA backbone is neutral be- 
cause each phosphate group is bonded to a hydrogen. 
In contrast in aqueous solution, the overall system neu- 
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FIG. 11: (Color online) Histograms showing the number of 
times, N, a certain molecular energy level is located within 
the energy interval given by the position and bin width of the 
corresponding histogram bar. The purple lines in panel (a) 
represent the HOMO and LUMO levels of the DNA strand 
in the absence of the solvent. The dashed lines represent the 
averages taken over the histograms. Panels (b)-(d) show the 
distribution of the highest (lowest) occupied (un-occupied) 
molecular levels localized on individual Guanine (Cytosine) 
bases. The Guanine bases labelled Gl, G2 etc are arranged 
from left to right with Guanine Gl at the 5^ end and G6 at 
the 3^ end. 



trality is ensured by including two Na+ counter-ions per 
base-pair. In order to describe the water solvation shell 
we perform empirical-potential MD simulations by using 
the NAMD2 package [58], parametrized with the standard 
additive CHARMM27 force field for nucleic acids [59 . 
For water we adopt the TIP3P model [60 that is modi- 
fied in the CHARMM force field and includes Lennard- 
Jones parameters for the hydrogen atoms as well as for 
oxygen [61 . The interaction between the H20 molecules 
and gold is treated as a 12-6 Lennard- Jones potential as 
implemented in CHARMM force field, with parameters 
for Au (e=0.039 kcal/mol and cr=2.934 A) taken from 
the literature [62 . Periodic boundary conditions are ap- 
plied with a cutoff of 12 Afor long-range interactions. 
The initial configuration of water molecules surrounding 
the pGpC 6 base pair DNA between the gold slabs is 
obtained using the VMD [63 solvation procedure. In or- 
der to maintain the size and shape of the cell constant. 




E-Ep(eV) 



FIG. 12: (Color online) Projected DOS for different compo- 
nents of the DNA and solvent system. The bold dark lines 
show the PDOS averaged over the 100 MD snapshots. PDOS 
from all the individual snap-shots are shown in gray. 



we perform simulations in the micro-canonical ensemble, 
with re-initialized velocities to 300 K for every 1000 time- 
steps, each time- step 2 fs. 

Since our main objective is to examine the effects of sol- 
vation on the underlying electronic structure of A-DNA 
helices, we fix the atomic positions of the DNA strands 
in order to avoid conformational rearrangements. Thus 
only the water molecules and the Na+ counter-ions are 
allowed to move during the MD simulation. After the 
initial equilibration of 1 ns we continue to run simulation 
for a total time of 20 ns. From the final configuration we 
extract the water solvation shell, the thickness of which 
is chosen so that both the back bone of the DNA helix as 
well as the counter-ions are adequately hydrated by wa- 
ter. In total, 211 water molecules and 12 Na+ counter- 
ions are included in subsequent simulation at 300 K. The 
trajectory is recorded every 4 ps from the initial equili- 
bration of 1 ns to a total simulation time of 20 ns. 

We find that at equilibrium, the hydrated Na+ 
counter- ions prefer to detach from the PO4 groups of the 
DNA and remain suspended among the water molecules 
while the DNA back-bone is also strongly hydrated. Thus 
the interaction between the counter-ions and the DNA 
bases is strongly screened in aqueous solution. 

In order to extract a statistical picture of the energy 
level dynamics induced by molecule-solvent interactions. 
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FIG. 13: (Color online) Zero-bias transmission coefficient and 
corresponding DOS for one representative time snapshot of 
the DNA and solvent system taken from classical MD simu- 
lations. 



we select 100 uniformly spaced snapshots of the MD 
trajectory and calculate the electronic structure of each 
snapshot. The spacing between consecutive snapshots is 
large enough (120 ps) that the correlation between snap- 
shots maybe assumed to be weak. In figure [To) which 
shows the planar average of electrostatic potential at zero 
bias plotted along the transport direction for all the snap 
shots considered, large fluctuations in the local electro- 
static potential across different shapshots are apparent. 
These fluctuations are caused by the dynamics of the po- 
lar solvent molecules and act to locally gate the electronic 
energy levels of the bases along the DNA strand. As a 
result, the energy eigenvalues corresponding to molecular 
orbitals localized on the bases spread across a range of 
energies. 



In figure 11 we show histograms representing the num- 
ber of times a certain molecular level is located in an 
energy range given by the bin width of the correspond- 
ing histogram. In figure 11 'a), the histograms for the 



HOMO and LUMO energies of the DNA strand are 
shown. Firstly we observe that the average HOMO level 
in the presence of the solvent is shifted down by approx- 
imately 1.5 eV relative to the HOMO of the molecule 
in the dry. The position of the average LUMO how- 
ever changes little in comparison to the dry case. This is 
similar to what found for benzene molecules attached to 
gold via thiol groups [56 . Furthermore, the HOMO and 
LUMO energies fluctuate about the mean value with a 
standard deviation of 0.36 eV and 0.28 eV respectively. It 
is worth noting that although the peak of the histogram 
for the HOMO level in the presence of the solvent is 
~1.5 eV lower than the HOMO level in vacuum, energies 
near the right edge of the histogram are relevant for the 
observed I-V characteristics. The energy at which the av- 
erage DOS for the occupied Guanine levels (see figure 
starts to be substantial is seen to be ~0.5 eV lower i. 
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the HOMO level in the dry. Furthermore, the HOMO 
levels discussed here correspond to neutral DNA strands 
and as such should correspond to vertical ionization po- 
tentials. If one also takes into account the possibility of 
solvent reorganization around the DNA [36^ ^ in the 
case of a hole charge being present on the strand, the ef- 
fective ionization potential can be lower. Also shown in 



figure 11 are the histograms corresponding to energy lev- 
els localized on specific Guanine and Cytosine bases. The 
large amplitude of the fluctuations relative to the energy 
spacing between the molecular levels in the dry means 
that unlike in the static case, states localized on different 
bases can cross each other. In fact while in the dry the 
HOMO was always localized on the Guanine base at the 
5' end of the DNA strand, the dynamic energy disorder 
induced by the solvent rearranges the order of the states 
in the HOMO manifold. This means that at a given mo- 
ment in time (i.e. for a given water configuration) any 
of the energy levels associated to the six Guanine bases 
can be the HOMO. There is even the possibility that the 
instantaneous HOMO may be delocalized across several 
Guanine bases whose energy levels are close to resonant. 
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TABLE I: Table showing the number of time-snapshots in 
which the HO energy level localized on the Guanine base from 
the row entry is higher than the HO energy level localized on 
the Guanine base in the column entry. 

In table [l] we quantify the fluctuations between the 
states of the HOMO manifold. We list the number of 
snapshots in which an energy level localized on a given 
Guanine base is higher than all those localized on an- 
other bases. The Guanine bases are enumerated as Gl 
to G6 in sequential fashion starting from the 5' end of 
the strand. For example, it is seen that the energy level 
of the base Gl at the 5' end of the strand is higher than 
the energy level on the base G2 which is next to Gl, only 
63 percent of the time. Likewise, the energy level cor- 
responding to the base G6 which is at the 3' end of the 
strand is higher than the level on the base Gl about 21 
percent of the time. From the the table it emerges that 
the energy levels associated to the terminal bases in the 
strand have a higher tendency of remaining at the edge of 
the HOMO manifold. Thus for instance the state corre- 
sponding to the Gl bases is more frequently the HOMO 
than any other level, and the state associated to the G6 
bases is more frequently the deeper level in the mani- 
fold. However, all the other energy levels fluctuate quite 
considerably. 

The effects of such fluctuations on the electronic struc- 



lan ture are discussed next. In figure 12 we show the PDOS, 
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averaged over the 100 snapshots considered, for different 
components of the system. In the same figure we dis- 
play ah the instantaneous DOS calculated at each time 
snapshot (gray lines merging in a shadow). We note that 
the net (average) effect of the solvent is to gate the occu- 
pied levels by shifting them to lower energies and effec- 
tively increasing the HOMO-LUMO gap. Note that the 
PDOS originating from the Na+ counterions is largely 
confined to the unoccupied levels and located higher in 
energy than the primary peaks in the unoccupied DOS 
projected onto the Cytosine and Guanine bases. Thus 
the Na+ counterions are not expected to play any sig- 
nificant role in either electron or hole transport through 
the DNA itself, with the exception of the fact that they 
contribute to the electrostatic disorder. 

In contrast, a large fraction of the PDOS for the sol- 
vent water molecules comes from the occupied levels and 
is peaked about 5 eV below the Fermi level. Importantly, 
however, a significant amount of the water DOS extends 
into the vicinity of the Guanine HOMO levels and there- 
fore water molecules might also contribute to the electron 
transport across DNA. In figure [13] we show the trans- 
mission coefficient as a function of energy and the cor- 
responding DOS for one time snapshot. The snapshot 
chosen is one in which the HO and LU molecular levels 
are very close to the average HOMO and LUMO val- 
ues. Qualitatively, there is very little difference between 
the transmission coefficients in the presence and in the 
absence of the solvent. However, we observe that the 
presence of the water molecules enhances the tunneling 
transmission probability around zero bias compared to 
that of the dry situation which suggests that the decay 
of the wave-function at energies in the DNA bandgap is 
slower in water than it is in vacuum. 



VII. DISCUSSION AND CONCLUSIONS 

The present work establishes a number of clear facts. 
Firstly, it shows that coherent tunneling across DNA 
strands attached to Au electrodes with realistic bond- 
ing geometries is through the base-pair stack and not 
across the DNA backbone. This, however, is far from be- 
ing band-like transport as speculated a long time ago on 
the basis of an analogy with aromatic salts [64 . In fact, 
although an idealized homogeneous and periodic DNA 
chain possesses narrow conduction and valence bands, 
finite DNA fragments are dipolar molecules and such 
bands fragment into a manifold of extremely localized 
states that extend over one or two neighbouring bases. 
Furthermore, intra base-pair interactions between com- 
plementary bases are weak with the result that molecular 
states can be classified almost exclusively as belonging to 
one base or the other. In pGpC DNA strands in partic- 
ular, the HOMO manifold has Guanine character while 
the LUMO has a Cytosine one. The coherent transport is 
then resonant tunneling-like across such localized states. 
Interestingly the transmission function is affected by the 



linker groups which attach the strand to the electrodes, 
even if those are connected to the backbone. 

We then find that the I-V characteristics are gapped 
with a conductance gap of the order of 2 Volt, in good 
agreement with controlled STS experiments. The trans- 
port is dominated by HOMO conduction and this essen- 
tially refiects the HOMO-LUMO gap and the alignment 
of the DNA molecular levels with the Fermi energy of 
the gold electrodes. Such an alignment has been carefully 
evaluated by calculating the IP and EA for the molecules 
with an accurate self-consistent DFT method and then 
by exporting the result to the transport calculations via 
the ASIC scheme. The agreement with experiments is 
however broken when one looks at the magnitude of the 
current at the onset voltage(~ 1 Volt). The current is 
several orders of magnitude smaller than what measured 
in experiments even for strands shorter than those mea- 
sured. This indicates that coherent tunneling is not the 
main mechanism for longitudinal transport in DNA even 
for short DNA fragments attached to electrodes. 

The picture changes drastically when one considers the 
effects of a few water solvation layers on the transport. 
Water molecules create a dynamical dipolar field that re- 
arranges the relative energy positions of the individual 
states of both the HOMO and LUMO manifolds. The 
average HOMO manifold is shifted down in energy by 
about 1 eV while the LUMO remains roughly in the same 
position. In this situation then both HOMO and LUMO 
transport may become relevant although the position of 
the DFT-LUMO is too low within our current approach 
and so a conclusive statement on bipolar transport is not 
possible based on our results. In addition we find that 
the presence of the water broadens further the HOMO 
and LUMO manifold, but more crucially it changes the 
position of the energy levels associated to any given in- 
dividual base by as much as 1 eV. 

Based on all these ab initio results we can formulate a 
mechanism for electron transfer in DNA. This occurs by 
means of incoherent resonant tunneling through narrow 
energy levels spatially localized at the base-pairs. Levels 
are brought into resonant by the dynamical water dipolar 
field, which then acts as a main source of de-phasing. We 
note that the typical time-scale for the re-arrangement of 
the DNA energy level due to the solvent is of the order 
of a few femto-seconds i.e. it is considerably faster than 
the DNA molecular vibration characteristic time. We 
then conclude that phonons play only a secondary role 
in DNA conduction, which in turn is determined by the 
dynamics of the solvent. Our calculations thus pave the 
way for the construction of quantitative models for DNA 
conduction based on a first principles Hamiltonian. 
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